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COMPUTERIZED TOMOGRAPHY USING VIDEO 
RECORDED FLUOROSCOPIC IMAGES 

A. C. Kak , C. V. Jakowatz, Jr,", N, A# Baily*^ and R, A, Keller^ 

In this paper we Investigate the possibnity of constructing computerized 
tomograms using data collected from a fluoroscopic system, it is shown that 
through proper handling of this data, useful images can be obtained. The sys- 
tem offers the advantages of eliminating the need for a highly stabilized, 
linear translating mechanism, and also of requiring relatively low patient 
dosage. In addition, the data gathering can be done in essentially real time. 

Several significant problems arise when fluoroscopic data, gathered via 
an Image intensif ier/vidicon combination, is used for tomographic purposes. 

First, it turns out that it is difficult in such a system to accurately determine 
a value for the reference x-ray intensity. This is due to scatter and to the 
fact that optimal use of the quantizer range, is desired. We have shown theoret- 
cally that the effects of such an error are to introduce a rIng-like artifact 
into the reconstructed tomogram, and have shown how this problem may be elim- 
inated. Second, a very important consideration with the fluoroscopic data Is 
that no more than sixteen gray levels can be meaningfully recorded directly by 
the vidicon. This low s ignal-to-nolse ratio would clearly not be suitable for 
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reconstructions of tissue density differences on the order of one percent. 
However, since a large number of points, say 500, may be digitized from a given 
TV line, and since a frame may be recorded in only a few milliseconds time, 
spatial and temporal averaging may be performed to reduce considerably the noise 
variance, without suffering significant loss of spatial resolution. This ail lows 
accurate reconstructions to be done. Finally, the fact that a diverging radi- 
ation source is employed must be taken into account. We have demonstrated that 
one can avoid having to either rearrange data or develop efficient algorithms 
that directly handle a fan beam geometry by simply using a parallel ray con- 
volution/back-projection algorithm, when the beam divergence is less than about 
15 °. Through computer simulations we have demonstrated these effects, in add- 
ition, we have reported results from our studies done on lucite phantoms and a 
freshly sacrificed rat. 
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1. INTRODUCTION 

The problem of digitally reconstructing the Internal structures of an 
object from measurements of Its two dimensional projections, resulting from 
transmission of radiation through the object, has been of interest In mathe- 
matics, radio-astronomy, and biological sciences for over 50 years. About 
four years ago this problem suddenly assumed major Importance In medical engi- 
neering after Hounsfleld and his associates [23] announced the development of 
a sophisticated computerized x-ray scanning system capable of reconstructing 
cross-sectional Images (tomograms) of the human head with high tissue-density 
discrimination capability. This development was heralded as a major break- 
through In diagnostic radiology because It, for the first time, permitted non- 
Invaslve visualization of the ventricles of the brain, as well as a large 
class of brain tumors and Injuries. Hounsfleld's scanning system, also known 
as the EMI scanner, does have one disadvantage, however. Each projection re- 
quired for constructing the tomogram Is obtained by linearly moving a highly 
collimated x-ray beam across the body section (Fig. 1). This means that 
movements of x-ray tube and the crystal detector have to be highly stabilized. 
This adds considerably to the cost of the system, which at present Is beyond 
the reach of many medical Institutions, This may also lead to poor quality 
tomograms If the object or parts of It are In excessive motion during each 
linear scan, such as may be the case In attempting to Image the heart. 

In this paper we have examined a computerized tomographic Imaging sys- 
tem which employs video- recorded fluoroscopic Images as Input data. By hook- 
ing the video recorder to a digital computer through a suitable Interface, 
such a system permits very rapid construction of tomograms. In fact, the 
system can be made In essence "real time”. With the use of fluoroscopic 
Images, the tomographic Imaging system Is essentially equivalent to the one 



shown in Fig. 2. Note that an entire projection is now recorded in a single 
Instant of time, reducing the total time required for all of the data collec- 
tion. Also, if the object now changes in time and if these changes are peri- 
odic, different projections can be taken at time instances synchronized with 
periodic variations, so that the reconstructed tomogram will represent the 
cross-sectional slice of the object at one instant of time. 

While computerized tomography using fluoroscopic data does have the ad- 
vantages of requiring lesser data collection time and being Immune to periodic 
variations in the object structure than the more familiar system requiring the 
linear motion, it does suffer from the disadvantage of the inherently low 
signal to noise ratio present in fluoroscopic images. As was shown earlier [ 3 ] 
the signal to noise ratio problem can be alleviated by suitable spatial and 
temporal averaging of the data. 

It is clear that both techniques for data collection, the linear motion 
of the source-detector in the EMI machine and the video recording of fluoro- 
scopic images, are equally important for computerized tomography. While the 
superior signal-to-noise ratio in the former are ideal for imaging the human 
head where good tissue discrimination capability is required, the latter tech- 
nique is superior for imaging parts like the thorax where periodic variations 
in the internal structure occur. 

In addition to the low signal to noise ratio, the other major distortions 
present In the fluoroscopic data arise from: the approximately logarithmic 

non-linearity Introduced by the fluoroscopic screen; errors in recording the 
air-transmission values, which result in the presence of ring-like artifacts 
when convolution algorithms are used for reconstruction; and the fact that 
projections are taken by non-paraliel radiation. A recent paper [3^] dealt 
with the problem of non-parallel radiation by using an algebraic reconstruction 



algorithm In which simultaneous equations were formed to explicitly take into 
account the non-parallel nature of rays. However, the algebraic techniques 
in general are known to be less accurate than the convolution techniques. 
Although the development of convolutional techniques for projections taken 
with radially diverging radiation has recently been reported , such algo- 
rithms seem to require large computer time. In this paper we have shown 
that for beams with divergence up to approximately 15 degrees one can use the 
convolution algorithm designed for the parallel radiation case with negligible 
degradation both quantitatively and from the point of view of visual quality, 
in fact, it may be shown that for very small angles of beam divergence a con- 
volutional algorithm designed for the parallel beam case gives superior re- 
sults and is computationally more efficient than either the algebraic techniques 
or the convolutional algorithms for radially diverging data. 

In the next section we will first very briefly review the' 1 Iterature 
dealing with digital reconstruction of objects from their projections and then 
discuss the digital implementation of the convolution algorithm. 


This statement must be taken with some caution. The speed and the accuracy* 
of an algorithm is dependent upon the pattern that needs to be reconstructed. 

Our conclusions are based on the comparisons made by many researchers on phan- 
toms relevant to x-ray computerized tomography In medicine. For example, 
for the phantoms of the human head section, Shepp and Logan l33] compared the 
parallel ray convolution algorithm with the algebraic technique and concluded 
that the algebraic algorithm required considerably more computing time to obtain 
a reconstruction of comparable accuracy. 
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2. RF.CQNSTRUCTION OF OBJECTS FROM THEIR PROJECTIONS 

Let f(x,y) be a two dimensional function shown in Fig. 3. 

Let the integration of the function f{x,y) along lines subtending 

angles of 6 + 90° with the x-axis be denoted by P^Ct) where t is a 

o 

distance parameter as shown in the figure. The P0(t) is called the 

projection of the function at an angle 6. The problem of 

reconstructing an object from its projections is then equivalent 

to reconstructing f(x,y) from a knowledge of p0(t) various 6. 

Many investigators Pf -9, 1 1-22, 24-30, 33»35] have addressed 

themselves to the above problem. A number of different approaches 

have been developed, the more important falling in three categories 

the algebraic methods [16-18,20-22], the Fourier methods, [11-15, 

26,28,35] and the convolution or the filtered back-projection 

methods [4,5,7“9»24,29,33] . Of these techniques, the convolution 

or the filtered back-projection methods are known to be computation 

ally most efficient and accurate. We will now discuss the digital 

implementation of the techniques that fall in this category. 

Let S„(u) be the Fourier transform of the projection Po{t), 

- 0 0 

that is 

Sg(u) - Og(t) (I) 

Given S_{u) for all 6, it is known that the two dimensional 
function f(x,y) (Fig. 3) can be reconstructed by using the follow- 
ing relationship 

f{x,y) = /q Q0(x Cos 6 + y Sin 9 )d 6 


( 2 ) 



7 


where 


Q„{t) = 03 (a) |u|e-i2™t 


du 


(3) 


The above formulas for reconstruction say that from each 
projection f* 0 (t) we calculate a "filtered projection" Q 0 (t) by 
using (1) and { 3 ), and use (2) to reconstruct the function f(x,y). 

The parameter u has the dimension of spatial frequency. The 
Integration In ( 3 ) must, In principle, be carried out over all the 
spatial frequencies. In practice the energy contained In the 
Fourier transform components above a certain frequency would be 
negligible. So for all practical purposes the projections may be 
considered to be bandllmited, an argument used by Mersereau and 
Oppenhelm [27] to arrive at some Interesting theoretical results. 
Next, let W be the smallest frequency beyond which the spectral 
energy In all the projections may be ignored. Then by using the 
sampling theorem [32], the projection can be represented by 




T p (JLn 

, ^ 0^2vr 

k=-<M 


Sin 2-nW(t~) 
2TrW(t-~) 


{^) 


Substituting (4) In (1), we get 


Sg{u) 


_L 

2W 


2 

k=-oo 




( 5 ) 


where 

by(u) = 1 hi 5. W 

= 0 


otherwi se. 
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Now if the projection functions are of finite order [27], which means 
they can be represented by N+1 samples for some value of N, ( 5 ) 
reduces to 

N/2 ■ k 

= 6 <“> - -k 


We will arbitrarily assume N to be an even number. Now the func' 
tions Sq(u) are zero outside the interval (-W,W) of the spatial 
frequency axis. Suppose that in this interval we desire to know 

i 

each'Sg(u) at a set of equispaced points gi-ven by 


u 


2W , 

m for m 


N 

2 ’ 


0 , 


. . . . , 


H 

2 * 


Substituting these in (6), we get 


c r 2Wx J_ p ,_k_v 


_ N n /7\ 

m "* 2^^ ••••} Uj m m m • j ^ V// 

The above equation is the familiar discrete Fourier transform (DFT) 
relationship and, therefore, can be rapidly evaluated by the fast 
Fourier transform (FFT) algorithms. 

Given the samples of a projection, the equation (7) gives the 
samples of Its Fourier transform. . The next step is to evaluate the 
"modified projection" Qci{t) digitally. Since the Fourier transforms 
S„(u) are bandlimited, (3) can be approximated by 

O 

Qg(t) = Sg(u) |u|e^^'^“^du 


(8) 
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N 


N/2 

S 

m=“N/2 


.. f 2Wv 


\rM\ 


J27rm~t 
e N 


(9) 


provided N Is large enough® Again, if we say that we would like 
to know the modified projections Q0(t) for only those t at which 
the projections P0(t) are sampled, we get 


-0 - (' 


2Wx 


N/2 

E 

m=-N/2 


S,^(m^) |m-rrl 


^1 

N 


mk 

eJ2-rrm^ 


k = “N/2, ,®o, “1, 0, 1, N/2 (10) 


By the above equation the function O^Ct) at the sampling points of the 
projection function is given by the inverse DFT of the product of 

O W 0 w 

S-(m~p) and im— |. By the familiar convolution theorem for the case 
of discrete transforms, (lO) can be written as 

^0 “ IT ’’e ^2V7^ * ^ ^ 


where * denotes the operation of circular convolution and where 


K I I ix 

^ IS the inverse DFT of the discrete function ^ - "r'j 


N 


'2W 


“10 1 — 

,, I, u, i, .o®, 2 “ 


Clearly at the sampling points of the projections, the functions 
Qg(t) may be obtained either in the Fourier domain by using (10), or 
in the signal domain fay using- (H),® The reconstructed picture f(x,y) 
may then be obtained by the discrete approximation to the integral in 
(2) , i ®e® , 


4 * 

In practice superior results are usually obtained by padding the pro- 
jection data with a sufficient number of zeroes so that for k = -N/2, 
CO®, N/2 the modified projection is equal to the aperiodic convolution 
of the projection data and filter impulse response. This is because 
in practice the assumptions of finite bandwidth and finite order are 
not strictly satisfied. 
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IT 

fU»y) =7 S. Q- (x Cos e. + y Sin e.) (12) 

where the K angles 6. are those for which the projections Pg(t) are 
known . 

Equation (12) calls for each filtered projection Qg (t) to be 
"back-projected”. This can be explained as follows. To every point 
(x,y) in the image plane there corresponds a value of t (= x Cos 9 
+ y Sin 6) for a given value of 6, The contribution that Qg, makes 
to the reconstruction at (x,y) is the value of for the corresponding 
value of to This is further illustrated in Fig, 4, It is easily shown 
that for the indicated angle 9., the value of t = (x Cos 9.+y Sin 9.) 
is same for all (x,y) on the line LM, Therefore, the filtered projec- 
tion will make the same contribution to the reconstruction at all 
these points. From this follows that in reconstruction each function 

Q„ (t) is smeared back over the image plane. The sum (multiplied by 
®i 

tt/K) of all such smearings results In the reconstructed image. 

Note that the value of x Cos 9. + y'Sin 9. in (12) may not cor- 
respond to one of the values of t for which Qg is determined in (10) 

°i 

or in (11), However, Q. for such t may be approximated by interpo- 

O • 

I 

lation. We have used linear interpolation, which naturally introduces 

some errors since Q. (t), in general', is not linearly dependent upon t, 

®i 

The reconstruction technqiie discussed here was tested by computer 
simulation on a number of test patterns. One such pattern is shown in 
Fig. 5a which is obtained by superimposing elli'pses. The advantage 
of using ellipses for computer simulation was first pointed out by 
Shepp and Logan [33]. This relies on the fact that one can write down 
simple analytic expressions for the projections of such patterns. 
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Note that projections obtained in this way are the "real “projections” 
as opposed to "pseudo-projections” [21 ]« In Figo 5a there are fine 
gray level differences between various regions. However, these cannot 
be seen in the photograph. Fig. 5c illustrates some of these differ- 
ences, The heavy dotted line in Fig. 5c shows gray levels across the 
middle horizontal section In Fig. 5a. 

Sixty projections of the patterns In Fig. 5a were calculated by 
the method in Appendix A, Each projection was sampled at 101 points. 
The Image reconstructed from this data by using equations (11) and (12) 
is shown in Fig, 5b, The light solid line In Fig. 5c shows the gray 
levels across the middle horizontal section In the reconstruction in 
Fig. 5b, 



12 


3. A SYSTEM FOR COMPUTERIZED TOMOGRAPHY USING FLOUROSCQPIC DATA 

Figure 6 shows a block schematic of the system used for the 
collection of the projection data. The phantom was mounted on a 
turntable. For each position of the turntable two frames of the 
projected image of the phantom were recorded on the video disc by 
the image-intensifier vldicon combination. Because of the high 
rate at which the data is produced by the vidlcon, it cannot be 
directly digitized. The use of the video disc and the scan con- 
verter is to slow down the data rate before its digitization by the 
A/D converter. 

Only a small area of each TV frame (Fig. 7) recorded on the 
video disc is actually digitized. The vertical width of this 
digitized area is approximately equal to the thickness of the cross- 
sectional slice of the phantom for which the tomogram is con- 
structed. As shown in Fig. 7» a strip of aluminum is attached to 
the face of the flouroscopic screen of the image-intensifier so 
that a part of the digitized data for each projection Is a measure 
of the x-ray transmission through aluminum. These values are used 
for scan-to-scan normalization of the data against variations In 
the x-ray generation. Note that the aluminum strip Is outside the 
image of the phantom on the image Intensifier. Aluminum Is used 
to prevent the saturation of the flourescent screen by x-rays that 
travel directly from the x-ray source to the screen. 

The part of each TV frame that is digitized, as shown in Fig. 7» 
is sampled on a 240 x 8 grid of points, 240 horizontally and 8 
vertically. After digitization, spatial averaging is performed over 
blocks of 3 X 8 points (3 horizontally and 8 vertically), so that 
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each 2^0 x 8 record reduces to a string of 80 data values. The 
regional averaging performed reduces such noise as that generated 
by quantum mottle, the image amplifier - TV chain, recording, 
playback and digitization. Temporal averaging further reduces the 
effects of nearly all types of noise. Temporal averaging was 
employed to a limited extent by recording two frames for each 
projection, and then averaging the results. 

Each of the 80 data elements for each projection was repre- 
sented by a 7 bit word, in other words, the quantization of the 
spatially and temporally averaged data was done to 128 levels. It 
should be noted that such a fine quantization has no meaning without 
the spatial and temporal averaging employed. That is because even 
under the best of laboratory conditions the flouroscopic data as 
directly recorded by the vidicon do not have more than fourteen 
gray levels [1 ,2, .10], implying a high statistical variance of the 
individual data points. As we have indicated later, the averaging 
employed reduces the noise variance of the data by a factor approx- 
imately equal to 48. - - 

By rotating the platform in Fig. 6 one- can take projections of 
the test object at the desired number of angles, usually between 60 
and 180. Of course, in a clinical version of the system, one would 
rotate the x-ray source and the detector combination around the 
patient, as opposed to rotating the object we have done. The two 
types of motions for generating projections are equivalent, however. 

In the next few sections we will by theory and/or computer 
simulation analyze the effects of some of the problems associated 
with flouroscopic data in the context of computerized tomography. 
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4. LOGRITHMIC NON-U NEAR IT I ES tN DATA AND THE EFFECTS OF ERRORS (N 
MEASUREMENT OF AIR-TRANSHISGION VALUES 

Consider a parallel beam of x-rays shown In Fig, 8a, Let be the In- 
tensity of the beam before transmission through the object. Evidently, Is 
proportional to the number of photons per unit area per unit time. Let us 
now consider the propagation of x-radlatlon along a particular ray labeled CD 
in Fig, 8a*, In the absence of any absorption the intensity at D Is also equal 
to Nq. However, in the presence of absorption by an object In the path of the 
beam, and under the assumption that the x-rays are monochromatic, the Inten- 
sity at D is given by the Larabert-Beer relationship 


N = 


g“/^a(s)ds 


(13) 


where ct(s) is the composite attenuation coefficient of the object along the 
path CD. 

In a fluoroscopic imaging system the beam of x-radiation is fan shaped 
rather than parallel. Therefore, the x-ray intensity decreases not only as a 
result of absorption In-the material, but also due to beam divergence. The 
x-ray beam for our system can be approximated by that shown In Fig. 8b. Let 
the total number of photons emitted by the source 0 per unit time Into the 


In practice the x-ray beams are polychromatic. The spectrum of a polychromatic 
beam changes continually as It propagates through an inhomogeneous medium. An 
Interesting discussion on the difficulty that this causes in interpreting the 
numbers In a computerized tomogram can be found in [36], Under the condition 
that the detector integrates al 1 the transmitted photons, it has not yet been pos- 
sible to incorporate the Interaction of a polychromatic beam with an inhomogeneous 
medium In the theory of reconstruction of objects from their projections. It Is 
for this reason various researchers in the theory of reconstruction have ^continued 
to make the assumption of monochromaticity. Note that in discussing their re- 
sults for a phantom consisting of a number of 2.5 cm rods of commerical grade 
materials in a fixed-length water bath and 120 KeV operation, the authors in [ 36 J 
point out that the numbers produced by the EMI scanner give very closely the 
attenuation coefficients calculated for monochromatic x-rays of energy 73 KeV. 

The fact that the assumption of monochromaticity here does not give meaningless 
results Is olear from the corroboration of the equation (31) by the ring like 
artifact in the tomogram derived from the experimental data in Fig. 17b. 
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cone OPQRS be N^. Let be the x-ray Intensity at a distance r in terms of 
the number of photons per unit area per unit time. Clearly, under the assump- 
tion that x-ray emission within the cone has no angular dependence, the follow- 
ing result follows from the conservation of the total number of photons: 

= Nq 


Equivalently, 


r ^2 


r #e 


(14) 


Note that under the idealization that the radiation source can be considered 
to be a point source, should approach infinity as r lends to zero. 

Equation (14) gives the attenuation due to beam divergence of x-ray in- 
tensity along a ray like OM in Fig. 8b. If now an absorbing material is 
placed in the path .of the ray OM as In Fig. 8c the intensity at the point M 
will be g i ven by 


N =-^ 


^-/”a(s)ds 


(15) 


where again Oi(s) is the attenuation coefficient of the material along the ray 
OM. For practical purposes it is more convenient to express the intensity at 
point M in terms of the intensi.ty at the point A in Fig. 8c for a ray like OA 
which does not pass through the object. Intensity values at points like A in 
Fig. 8c are called air-transmission intensities. Let be the air-transmission 
Intensity. By (14), Nq = r^ . Substituting this In (15), we get for the 
intensity at M 

N - N. 
r A 


( 16 ) 
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The composite attenuation coefficient a Is In general a function of 
position In the absorbing material, that is, a function of the coordinates 
X, y, and z shown in Fig, 8c. If the angle 6 is small enough, the dependence 
on z can be ignored and the attenuation can be considered to be essentially a 
two dimensional function a(x,y). Expressing this dependence explicitly, (16) 
can be written as 

./• a(x,y)ds 

e Vay OM (17) 

Since In the system described in Section 3, the measurements in the strip 
PQRS in Fig. 8c are integrated along the shorter direction, we have essentially 
the two dimensional system shown in Fig, 8d, where we have ignored the curva- 
ture of the strip PQ,RS. It is clear from the above discussion that the mea- 
surement ! (t) at, say M”^, which is a distance t from the indicated reference 
point. O'*, is given by 

-/a(x,y)ds 

I (t) = i^^e L(t) (18) 

where 1^^ is the air-transmisslon measurement at a point like A"* and where 
L(t) Is the ray path from the source 0 to M”*. 

The above relationship is valid only In the absence of any non-linearities 
Introduced by the detection mechanism. However, it has been experimentally 
observed that the recorded measurements -a re approximately a logrithmic func- 
tion of the x-ray Intensity impinging on the Image intensifier tube (see Fig. 

* 

9). Therefore, over a large dynamic range what Is actually measured Is a 
quantity proportional to lln l(t). We get from (18) 

£n l(t) = -/j^^^ja(x,y)ds + £n I^^ (19) 

— _ 

Further discussion on this can be found In "Response of Image Intensifier 
Fluoroscopic - TV System," by Bally and Keller to appear In INVESTIGATIVE 
RADIOLOGY. 
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From FIga 8d it is clear that /j^^^ja(x,y)ds is the projection of the 
function a(x,y) by a diverging beam on the line Evidently, 

projections at different angles may be taken by rotating the source- 
detector arrangement around the object. In the preceding section we used 
Pg(t) to denote the projection at angle 6. From (19), we get 

pQ(t) = - An 1 (t) + An (20) 

where the data An l(t), of course, depends upon the angle 0 at which 
the projection is taken,. 

It is clear from (20) that if Pg(t) is to be known correctly, the 
air-transmission values An 1^^ must be known. Experience has shown 
that in practice 1^^, if taken simply from a ray like OA"* in Fig, 8c, 
is in general not known sufficiently accurately, the errors being large 
enough to cause the presence of structured artifacts in the tomograms. 

The errors in the measurement of may be caused by one or both 

of the following factors. 

(i) Limitations of the quantizer in the digitization process 

[32]: The output of!_the scan converter (see Section 3) consists of 

samples that may take any value in a continuous range. The quantizer 

divides the continuous range into a fixed number of intervals and 

all the samples with values in an interval are assigned a fixed 

value at the output of the quantizer (Fig, lO), Evidently, the 

operation of the quantizer requires that the lower and the upper 

limits, I . and I , respectively, of the continuous range be 

specified. To minimize the quantization error for a given number 

of output levels, I must be close to the largest value and 
’ max 

I . close to the smallest value of the Intensities that 
mi n 

result after transmission through the object. If I is 

max 
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chosen close to the air-transmission values, all the useful 
information (i.e.,rays that do go through the object) would get 
crowded toward the bottom end of the quantization scale, leading 
to large quantization distortion. This implies that the air- 
transmission values may in practice be much larger than I 

max 

Since for any value larger than I , the output of the quantizer 

max ^ 

is q|^(Fig,10) ^ there may be greater errors in recording air- 
transmission values than other data. 

ii) Scattering: The interation between photons and matter 

results in the removal of the photons from the original beam 
either by scattering or by absorption. The attenuation of an 
x-ray beam caused by both these phenomena is represented by the 
coefficient a in equations (13) through (19). Some, scattered 
photons reach the flourescent screen of the image intensifier 
and cause the degradation of the signal to noise ratio there. Now 
if to reduce the noise variance, the air-transmission values are 
calculated by averaging the data for rays that do not pass through 
the object (after the field inhomogeneities have been subtracted 
out), it is clear that due to the non-negative nature of the con- 
tribution by scattering, this average will be greater than that in 
the absence of scattering. 

We will now show theoretically that the errors in the air- 
transmission data (5-n l^.. in (20)) will lead to the presence of 
ring like structures in the reconstructed tomograms. 

The equation (3) for filtered projections, Q-(t), can be 


written as 
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Qg(t) - (u) [-j sgn(u)]e-^^™^du 

where the function sgn(u) denotes 
sgn (u) = 1 for u _> 0 
= -1 for u < 0 


( 21 ) 


( 22 ) 


By using the standard convolution theorem (21) can be written as 

Qg(t) ={i.F.T. of j2iru Sq(u)}*{I.F.T of sgn(u)} (23) 


where the symbol * denotes convolution* and the abbreviation I.F.T. 

0 

stands for inverse Fourier transform. The I.F.T. of j2TruSg(u) is ■^F0(^) 
while the I.F.T. of sgn(u) is y* 'Therefore, (23) can be written as 


3P.(t) 

- -It- 


2lT^t 


(24) 


Now, the second term on the R.H.S. in (20) is not a function 
of t. Also, note the air-transmission value, once it is determined. 


is same for all the projections. (The data is corrected for scan- 
to-scan variations.) Therefore, errors in £n cause a d.c. 
shift in all of the projections. Let this d.c. shift be denoted 


by d^. If Po(t) denotes the correct projection and P''p(t) the 

'a' 


e 6'"' vw- .... f-.-j . 

recorded projection in the presence of errors £n I^, it is clear 




(25) 


Since the reconstruction algorithm discussed in Section 2 is linear, 
the tomogram made from projections P'*Q(t) obtained from the recorded 
data will be the sum of tomograms from the Pg(t) and whatever 
results when each projection is taken to be a constant equal to d^. 
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Therefore, the distortion in the reconstructed tomogram is 
another tomographic image caused by d In (25). If e(x,y) denotes 
the distortions in the tomogram, ciearly from (2) and (2^f) 


e(x,y) = /J Cos 0 + y Sin 0)de 


( 26 ) 


where 




(27) 


27T t 


For the purpose of illustration let us assume that the object 
is bounded by a circle of unit radius and that the data for all the 
projections is recorded over a length of 2 (i.e., for t between -1 
and +1 ) . The d.c, shift in each projection is then a rectangular 
pulse as shown in Fig. 11a, The derivative of this function is shown 
in Fig. lib and can be expressed as 


de = d^[6(t+l) - 6{t-l)] 


( 28 ) 


Substituting this in-(27), we get 


Q-e(t) 


d 


/(1-t^)' 


(29) 


Substituting (29) in (26), we get for the distortion e(x,y) in 
the tomogram 


^ ^ 2d 

e(x,y) = — Y /J 2 

2ir l**(x Cos 0 + y Sin 0) 


d0 


(30) 


The integration in (30) can be simplified by noting that the 
distortion caused by the error in the air-transmission value should 
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be a circularly symmetric function. That is, if r is the radial 

distance from the center of the picture, e(x,y) should be function 

/2 ? 

of the form E(r) where r = A +y . This function E(r) can be 
obtained by putting x=y in (30) and writing the L.H.S, as E(r). 

We get: 

2d 


“ 2 -^0 , 2 


de 


2tt 


1 - X (Cos 0 + Sin e)' 


w I 


(31) 


The above result shows that the error in the air-transmission 
value should result in the appearance of circularly symmetric artifact 
in the tomogram, the pattern being such that its intensity increases 
sharply near the outer edges of the tomogram. Even though in theory 
these ring like artifacts take infinitive values at the outer 
boundary of the picture, in actual digital implementation these will 
be finite, their magnitude depending upon the spatial sampling density 
in the picture plane. 

In order to simulate the effects of errors in the air-transmission 
values, a constant value was added to all the projections of the pattern 
in Fig. 5a. Fig, 12 shows the reconstructed image when the value of 
this constant is approximately 30 percent of the maximum gray level 
difference in Fig', 5a. Note the circular artifact outside the main 
pattern in Fig, 12, This is In keeping with the theoretical conclusions 
just arrived at. Also note that even though from a qualitative (visual) 
standpoint the artifacts caused by air-transmission errors lie primarily 
outside the main pattern in Fig. 12, from the quantitative standpoint 
the distortion exists over the entire imagej its form given by eq. (31). 
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In section 7 artifacts of the type discussed here will be shown 
to occur in the tomograms of a phantom and of a thoracic section of a 
rat. Based on the discussion presented here, we will remove these 
artifacts by introducing a correction in the air~transmission value. 
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5. QUANTIZATION EFFECTS 

In general, the noise level present In the TV and video disc recorder, 
even under the best of laboratory conditions, allows a maximum of fourteen 
gray levels [1-3,10]. This represents a high noise content. If an attempt is 
made to achieve reconstructions from this data without any preprocessing, the 
resulting pictures would generally be of low quality. 

To get around this problem of low signal to noise ratio, the spatial and 
temporal averaging described in Section 3 was employed. Each data element in 
the projections is obtained by spatially averaging 2k samples of the recorded 
images and temporally averaging over two recordings. With this averaging, the 
noise variance of the projection data should be approximately 7^ of the orig- 
inal noise variance, the result being strictly correct only if we can assume 
that the noise in each image sample is independent and identically distributed. 

Because of the considerable reduction in the noise variance, the data 
obtained after spatial and temporal averaging should contain many more useful 
gray levels than the dozen or so in the fluoroscopic image. It was because of 
this reason that a 7“bit~quantizer was used, yielding 128 gray levels for the 
projection data. 

The discussion here does bring to light the fact that if fluoroscopic 
data is used for reconstruction, the degree of quantization is dependent upon 
the extent of the spatial and temporal averaging employed. 

It is clear that the degree of quantization has to be matched to the 
signal-to-noise ratio of the data. While over-quantization, in general, re- 
sults in increased demands on the core and other peripheral memory requirements 
and, sometimes, reduced speeds in A/D conversion, under-quantization always 
leads to deterioration of the reconstructed image. To illustrate this effect 
by computer simulation, the samples of the analytic functions for the 
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projections of the picture in Fig. 5a were quantized to 8 levels. The recon- 
structed tomogram with this quantization is shown in Fig. 13a. When the data 
is quantized to 16 levels, the reconstruction is shown in Fig. 13b. Note that 


this is approximately the picture that would result if in a^/fluoroscopic sys- 
tem no spatial or temporal averaging was employed becaus^-^e number of useful 
levels in the data would only be approximately 16. When the data is quantized 
to 32 and 128 levels, respectively, the reconstructed images are shown in Figs, 


13c and 13d. 


We have said approximately because while the real fluoroscopic data has random 
noise, the quantization noise in the simulated data is more structured. The 
structured noise may generate artifacts by interference. A more realistic model 
would be to add random noise to the simulated data before or after quantization. 
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6. EFFECTS OF DIVERGING RADIATION 

In the system discussed in Section 3 the projections are taken 
with diverging radiation. One can take two approaches in con- 
structing tomograms from such projections: (I) Rearrange the pro- 

jection data in such a way that the new data is approximately what 
would be generated by a parallel beam of radiation, and then use 
an algorithm like that In Section 2 to make the tomogram; or 
(ii) Develop an algorithm that directly yields a tomogram from 
the projection data taken with diverging radiation. 

In Fig. 14 we have shown how one may rearrange the projection 
data obtained with diverging radiation to generate approximately 
the projections that would result if a parallel beam source was 
used. For the position 0, of the source, the projection data is 
recorded on the line D^D^. A different projection is obtained when 
the source is moved to O' and the detection line to The ray 

OA in the first projection is parallel to the 0'*A'* in the second 
projection. By grouping from all the recorded projections the rays 
that are parallel to, say, OA in Fig. 14 one can approximately 
obtain projection data that would be generated by a parallel beam 
of x-rays, parallel to the OA line. One can then use the algorithm 
in Section 2 to construct the tomogram. The disadvantage of this 
technique is the additional computer time required for rearranging 
the data. 

Clearly, it is desirable to have an algorithm that would di- 
rectly convert the projections obtained with diverging radiation 
into tomograms. It is relatively easy to devise algebraic types of 
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algorithms for this purpose [31]. However, algebraic reconstruc- 
tions are inherently less accurate and computationally less effi- 
cient than the convolution type algorithms such as that discussed 
in Section 2. Convolution type algorithms for data obtained with 
diverging radiation are still under development at Purdue Univer- 
sity and elsewhere and are not yet available for experimental 
work. 

We will now show by computer simulation that for upto 15° in 
beam divergence, if the convolution algorithm in Section 2 is used 
for constructing tomograms from projection data obtained with 
diverging radiation, the distortion introduced is negligible. For 
the pattern shown in Fig. 5a the computer was programmed to generate 
projection data using diverging rays by the method discussed In 
Appendix A. The projection data was generated for the following 
cases of beam divergence: 5°, 10°, 15°, 20°, and 60^. In each 
case 60 projections were determined and each projection sampled at 
61 points. The tomograms were constructed on a 60 x 60 array of 
points by the method in Section 2. The results are shown in 
Figs. 1 5a-e. Comparison of Figs. 15a-d with Fig. 5b clearly 
demonstrates thal^ for the reconstruction algorithm used, beam 
divergence up to 15 ° can be ignored. 
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7. EXPERIMENTAL RESULTS 

A test phantom is shown In Fig, 16, It consists of a 2 cm x 
5 cm X 7 cm Incite block taped to a one millimeter aluminum sheet. 

There are two holes in the lucite block, each of diameter 1 cm. The 
holes were filled with a 1.5% (weight/unit volume) solution of 
Hypaque. The contrast between the solutions and the surrounding 
incite Is provided by the K absorption edge of iodine. The phan- 
tom was mounted on a protractor which revolved about its center. The 
center of the block was also the center of rotation. This assembly 
was then placed midway between the x-ray target and the input 
screen of the flouroscopic system. The tube was set at approx- 
imately 60 kVp and the effective energy determined from absorption 
in aluminum was 22.8 keV. The tube current was approximately 1 mA. 

The beam divergence angle was estimated to be 10°. Sixty projec- 
tions at 3° intervals were taken and each projection digitized and 
spatially and temporally averaged, as discussed in Section 3- 
X-ray field inhomogeities were removed from the data by first 
recording a blank projection (i.e., measuring the x-ray Intensities 
without any object between the source and the image intensifier) 
at low tube current to prevent the saturation, computing the 
angular dependence of radiation and then correcting the projection 
data accordingly. 

Figs. I7a-c show the tomograms obtained for this phantom. Fig. 

16a is obtained when the air-transmission value is set equal to the lar- 
gest value at points outside the geometric shadow cast by the test 
phantom on the flourescent screen of the image- intensifier tube. It 
is clear from the theory and the computer simulation presented in 
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Section 4 that the ring artifact in Fig. 17a is due to the errors 
in the air-transmission value. The correct value of the air-trans- 
mission intensity can only be arrived at by trial and error, 
keeping in mind the fact that the correct value should make the 
ring in Fig. 17a disappear. Fig. 17b is the tomogram obtained 
from the same data as that used for Fig. 17a, except that the air- 
transmission value used has been corrected. By gray level slicing a 
tomogram one can display separately regions of different densities. 

The above experiment with a slight change was repeated with the 
test phantom replaced by a dead rat. Instead of rotating the platform 
In Fig. 6, the x-ray tube and the image intensifier were rotated 
around the rat, thus representing more nearly the clinical situation. 

A total of 180 equally spaced projections were taken through part of 
the thorax containing the heart. The approximate location of the 
section for which the tomogram was made is shown in Fig. I8a. Fig. 18b 
shows the constructed tomogram with the air-transmission value set 
equal to the data obtained from rays outside the geometric shadow 
cast by the rat body. The ring artifact resulting from the errors 
in the air-transmission value is again evident. Fig. 18c results when 
the air-transmission value is corrected. Fig. l8d identifies the var- 
ious anatomical features for the tomogram. 
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APPENDIX A 


DERIVATION OF EXPRESSION FOR PROJECTION 
DATA USED IN COMPUTER SIMULATIONS 


In this appendix we show how an analytic expression for the 
projection data, p0(t), for a certain class of picture functions may 
be derived. This data can easily be generated and subsequently used 
in the reconstruction algorithm in order to test the various effects 
addressed in this paper. 

Now it should be noted here that it is fairly easy to obtain 
a sort of “pseudo-projection data" for any digital picture l2l]. 

In order to do so, one can think of laying down "strips" centered 
around each line L(t,9) in a given projection, as shown in Figure 
Al.l. Now if the picture to be reconstructed assumes a constant 


va 


lue f.. in each square centered at (i,j),one can estimate P0(t) 


U 


by summing f.j for those block centers (i,j) which are included 
in the proper strip. This type of data, however, does not 
necessarily closely approximate the actual P0(t) 3nd hence is not 
generally very useful for testing the reconstruction algorithms. 

Now one could easily argue that this procedure could be made exact 
by simply taking the actual fraction of each square included by a 
given ray into account, rather than by simply summing the f.^ for 
the block centers (i,j) included. Although this can be done, it is 
a time consuming procedure to compute the appropriate fractions for 
each block for each ray of each projection. What is desired is an 
analytic expression for P„(t) for some class of pictures, so that 
exact test data could be rapidly generated. 
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Now as it turns out, the line integrals P^Ct) can be worked 
out precisely without difficulty if the picture function f{x,y) 
is such that it assumes a constant values, say C^, inside of an 
ellipse and is zero elsewhere. This was suggested in a recent 
paper by Shepp and Logan [33 3 • It will then be easy to obtain 
the projection data for any picture which is a superposition of 
such ellipses. To make things more precise, consider first an 
ellipse centered at the origin with semi-major and semi-minor 
axes given by A and B, respectively: 

2 2 

A B 

Now define f(x,y) such that: 

2 2 

f (x,y) = C. for ~ + £ 1 

‘ A B 

f(x,y) = 0 , otherwise 

In order to obtain the projection function P 0 (t) for f(x,y), 
we need to compute: 

for each line L(t,0) in the 0-projection, as shown in Figure Ai.2. 

[Note also that we further arbitrarily restrict the size of the ellipse 
to be such that it is contained entirely within the unit circle 
centered at the origin. In this way, ali of our simulated pictures 
will be such that f(x,y) = 0 for x + y >1. This, in turn, 
requires that P 0 (t) be zero for all jt| > 1 , for all 0 .] 
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Since f(x,y) = for all points inside the ellipse, it is 
clear that for a given line L(t,0) the value of Pq(t) = 
is given simply by times the length, Q.^^» of the line segment, 
Q^Q .2 Figure A1.2). If (x^y^) and (x^,y 2 ) are the coordinates 

of points and Q^, respectively, then it may be seen from 
Figure A1 . 1 that 

___ 1^2"^ J 

^1^2 sin e (A2) 


So we need to compute the points of intersection, Q,^ , 
line L('t,6) with the ellipse. 

Now the line L(t, 6) is given by: 

X cos 0+ysin8=T (A3) 


Substituting the value of y as given by (A3) in (Al) and calculating 
the roots of the resulting equation, we get 


X , = -5 p . sine e.B^sin=^e-T^)] 

(A^ cos 0 + B sin*^ 0) 
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Now (A4) implies that the value of x which forces the radical to be 
zero must be the value of x for which ~ and hence for which 
“ x^ I = 0. This is precisely the situation for which L(x,6) 
is just tangent to the ellipse, and so this value of x is the 
"projection half-width". This is denoted by a = a(0), and is 
illustrated in Figure A1.2. Clearly, we have: 

a(0) = cos^ 0 + sin^ 0 (A5) 

So we now obtain: 




2 AB sin 0 


/i 1 

/a “ t 


and substituting into (A2) yields: 


2C, AB 


/I T 

/a - t 


(A6) 


(A7) 


As a special check of the result ( A7 ) , it is clear that 
PQ{t) should assume a value of 2C^B for the case 0=0, t=0. This is 
evident from Figure A1.2. 

In order to obtain P„(t) when the ellipse described above is 
centered at (h,k) and is rotated through angle a , as depicted 
in Figure A1.3» only a slight modification is needed. 

It is not difficult to show that the projection data for this 
new picture function is related to P(t,0) of (A7) by: 


P**(t,0) = P[t - s cos (y ~ 0) > 9 ■ a] 


(A8) 


where. 


s = /h^ + k^ and y = ^ (•^) 
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That is, (s,y) are the polar coordinates of the translated center, 
(h,k). Finally, if a picture function f^Cx.y) is formed by 
superimposing N single ellipse pictures, where each ellipse assumes 
some constant value C. inside and is zero outside, it is clear 
that (All) may be modified to write the projection data for f 2 (x,y) 
as: 

N 

Pf (t,e) - 2 C. • P. [t - s. cos (y, •• e), 0 - a.) (A9) 

2 1=1 ‘ ' ' 

where s., y., and a. describe the position and orientation of each 

separate ellipse in the picture plane. 

For the case of fan beam (i.e., diverging beam) geometry, the 
same procedure may be followed, except that now the distance to be 
calculated in each case is as shown in Figure Al.A. 

Finally, it should be noted that no x-ray tube in practice has 
a zero spot size and no beam, however, collimated, will, have zero 
"thickness". It is, therefore, desireable to be able to simulate 
data that would result from "rays" of some thickness, as opposed to 
the ideal "pencil beams" considered thus far. To this end, consider 
the situation shown in Figure A1.5. Here it is assumed that the 
detectors measure the average value of the line integral of f(x,y) 
over the range of values, t^ - T/2 < t < t^ + T/2, where T is the 
"strip width". This, of course, is equivalent to finding the area 
of a given strip and dividing i t by T. For example, the projection 
data for this strip integral case may be found simply from A7 for 
the single ellipse xVa^ + y^B^ = 1 as: 
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P 


strip 



.t+T/2 

t“T/2 


2 ABC 

P(X,0) dX = =- 

a^T 


,t+T/2 

•'t-T/2 



dX 


AB C, 


2t 

a T 


[X»4^ - X^ + 


. -1 

s I n 



t+T/2 

t-T/2 


(AlO) 


Now the projection data for the case of N ellipses with various 
positions and orientations can again be found by analogy to (A9) . 
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FIGURE CAPTIONS 

1. A schematic diagram of the linear motion of the radiation source- 
detector combination required to generate each projection. 

2, The schematic diagram shows that in a flouroscopic system an 
entire projection is taken at the same time. 


3. The one dimensional function Pa(t) is the projection of the function 

in the direction of 0 + 90°. 

4. This figure illustrates the back-projection of the filtered pro- 
jection function Q,g(t). 

5. a) The pattern used for computer simulation. 

b) The image reconstructed from the projections of Fig. 5a. 

c) The relative gray levels in Fig. 5a. 

d) The relative gray levels in Fig. 5b. 

6. A block schematic of a system for computerized tomography imaging 
using flouroscopic images. 

7. The figure shows the part of the TV frame that is actually digitized. 

8. a) Absorbing material in a parallel beam of x-rays. This would 

be the case for linear-scan tomography. 

b) x-ray beam in a flouroscopic system is diverging as shown here. 

c) Absorbing material in a diverging beam. 

d) For small angles 0 in Fig. 8c the dependence on z-axis can be 

ignored and the resulting geometry is as shown here. 

9. This figure shows that the digitized data is proportional to log 1 
rather than 1, where I is the intensity of x-rays impinging on the 
flourescent screen. In our experiments we used a 7“bit quantizer 
(128 levels) and the projection data was typically contained in the 
range 30 to 120. For this range of the quantizer the curve is 

1 i near. 

10. All the analog values between z. and z._j^^ are assigned the discrete 
value q. by the quantizer. 

11. a) The d.c. shift In all the projections caused by the error in the 

air-transmission value after the projection have been corrected 
for scan-to-scan variations. 

b) Derivative of the function is Fig. 11a. 
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12. Computer simulation of the effect of air-transmission value errors. 
Note the ring artifact. 

13» These reconstructions indirectly show the importance of spatial 
and temporal averaging in a flouroscopic system, since it doesn'.t 
pay to have more quantization levels than that decided by the 
signal-to-noise ratio in the projections. The figure shows the 
reconstructions with projection data quantized to .(a) 8 levels, 

(b) 16 levels, (c) 32 levels, and (d) 128 levels. 

14. Figure shows how the data from diverging radiation can be rearranged 
to approximately look like the data from parallel radiation. 

15* These reconstructions show that the beam divergence effects are 

Insignificant for divergence angles up to 15°. The reconstructions 
shown are for (a) 5°, (b) 10°, (c) 15°, (d) 20°, and (e) 60°. 

16, (a) Lucite block with two holes that are filled with iodine 

solution. 

(b) The block is then taped to a 1 mm thick aluminum backing 
and mounted between two aluminum blocks. 

17, (a) Image of a tomogram through the lucite block. The ring is an 

artifact generated by errors in the air-transmission value. 

(b) The tomogram after the air-transmission value is corrected. 

18, (a) The approximate location of the thoracic section for which the 

tomogram was made. 

(b) The tomogram with no correction for the air-transmission value. 

(c) .Tomogram with the air-transmission value corrected. 

(d) Identification of feature in Fig. 18c. 








FIGURE 2 








Figure 5^ 






Gray Levels 



Fig. 5c 



A X-RAY SOURCE 

»«»**!» 

***•*♦, 




FIGURE 6 






ALUMINUM STRIP ON THE FACE OF 
THE IMAGE “INTENSIFIER TO PROVIDE 
SCAN-TO-SCAN NORMILIZATION 



FIGURE 7 





GRAY SCALE 
(QUANTIZER OUTPUT) 



2.6 


2.8 


FIGURE 10 


SAMPLE VALUE RANGE 



DECISION LEVEL 

OUTPUT LEVELS OF 
THE QUANTIZER 





FIGURE n 


Figure 12 





Figure 13b 







Flgur* 13d 








Figure 15b 










I 


59 



Figure 16b 




Figure I7a 





I 



FIGURE 18a 


62 














FIGURE A 1.3 



FIGURE A1.4 


^ ^ 3 



FIGURE A1.5 


